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1 . I nt roduct ion 



Finite difference methods specify the dependent variables at 



certain grid points in space and time, and the derivatives in the equations 
are evaluated using Taylor series approximations. The Galerkin procedure, 
which will be treated in this chapter, represents the dependent variables with 
a sum of functions which have a prescribed spatial structure. The coefficient 
associated with each function is normally a function of time. This procedure 
transforms a partial differential equation into a set of ordinary differential 
equations for the coefficients. These equations are usually solved with 
finite differences in time. The two most useful Galerkin methods are the 
spectral method and the finite element method. The spectral method, which 
employs orthogonal functions, has been used in meteorological problems for 
a number of years. The finite element method employs functions which are 
zero except in a limited region where they are low order polynomials. This 
method, which was developed in engineering, has only recently been introduced 
into meteorology and oceanography. 

The Galerkin procedure can be illustrated with the following equation: 



is a specified forcing function. Suppose that (l) is to be solved in the 
domain a _< x _< b and that appropriate boundary conditions are provided. Con- 
sider a series of linearly independent functions cp^ (x) which will be called 
basis functions. The next step is to expand u(x) into a series as follows: 



(u) = f (x) 



( 1 ) 




u is the dependent variable and f(x) 



N 



u(x) 




( 2 ) 
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where u . 

J 



is the coefficient for jth basis function. The error in satis- 



fying the differential equation (1) with the N terms of the sum (2) is 

N 

e N = X ( u^cp.) " ' ( 3 ) 

The Galerkin procedure requires that the error be orthogonal to each basis 
function in the following sense: 

/ ■ ■ " ■ * * ■ 



The final form is obtained by substituting (3) into (4): 



b N 



b 




cp^f(x)dx = 0 , i=l,...,N . 



(5) 



a a 

This reduces the problem to N algebraic equations which relate the unknown 
coefficients u^ to the M transf orms" of the forcing function. This procedure 
is quite general and can be applied to more dependent and independent 
variables . 

2 Example with Spectral and Finite Element Methods 

Now the spectral method and the finite element method will be applied 

to the following simple form of (l) : 

,2 

— ~ = f(x),o£x_<TT. (6) 

dx 

The boundary conditions are 



u(o) = u(tt) = 0 . 



(7) 



For the spectral method the following basis functions are appropriate: 



9 



epj = sin jx , j=l N . 



( 8 ) 



These functions are orthogonal on the interval 0 _< x _< tt and they satisfy 
the boundary conditions (7) . With these basis functions 

N N 

£<E 'v” ■ E ( - j2) v • 



j=l 



j-1 



and ( 5) becomes 
N 



TT TT 

.dx = / cp^f(x)dx , (9) 



- y vv* ’j 



The product of the basis functions can be written 
TT TT 

sinix sinjx dx = 



■*/ 



[cos(i-j)x-cos(i+j)x]dx = (tt/ 2) 6 , 



( 10 ) 



where 6.. is the Kronecker delta which satisfies 6^ . = 1 if i = i and 
ij ij 

S^. =0 if i 4 j • Equation (10) is merely the orthogonality condition 
which arises since the integral vanishes except when i = j . With the use 
of (10), the solution to (9) becomes 

TT 



u. 





cp^f dx . 



(id 



o 

Each coefficient is proportional to the finite Fourier transform of the forc- 
ing term. In this example both the error in the solution and the error in 
the differential equation are orthogonal to the basis functions. This is 
because is proportional to CfE so that if the error is orthogonal to 

cp ) it will also be orthogonal to cp. . This will also be true when certain 
l i 

other linear equations are treated with the spectral method, but it will not 

generally be true with nonlinear equations. 
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Now consider the same differential equation (6) with the finite ele- 
ment method. Divide the interval 0 £ x £ tt into N+l segments such that 
(N+l) Ax = TT . The basis functions are chosen to be tent shaped piecewise 
linear functions which are also called chapeau functions, as shown in 
Figure 1. As can be seen from the 



l 4 







(j-l)Ax jAx 

x -* 



(j+1) Ax 

Fig. 1. Piecewise linear basis function. 



ir=(N+l) Ax 



Figure, cfE (x) has a maximum of 1 at x = jAx , which is called the nodal 
point . The basis function decreases linearly to zero at x = (j-l)Ax , and 
it is zero everywhere else. Mathematically cp^ (x) is defined as follows: 

0 , x > (j+l)Ax or x < (j-l)Ax 
cp.(x) = j (x-(j-l)Ax)/Ax , (j-l)Ax £ x £ jAx] 

( (jfl)Ax-x) /Ax , jAx £ x £ (j+1) Ax 



( 12 ) 



Note that the coefficient u^ is actually the value of the function at 
x = jAx since cp. (jAx) = 1 and cp^(jAx) = 0 for i ^ j . These elements 
are not quite orthogonal, but only adjacent elements interact. The boundary 
conditions (7) are automatically satisfied although this is not necessary 
in many cases with finite elements. 

Equation (5) now becomes 



E 

4 = 1 




CPi 



f (x)dx = 0 . 



This form of the equation is not appropriate because it involves a second 

derivative of the basis function which is only piecewise linear. However, 

this problem can be avoided by integrating the first term by parts as follows 
7T 77 




d C Pi dc Pi dc Pi dc Pi 

L dx^ dx dx dx 



] dx -I 



CP± 



f(x)dx = 0 



The first term vanishes because all of the cp's. 

The Galerkin equation now becomes 
TT 77 




dcp. 
dx dx 



dx 




cp i f (x)dx 



i-1,. 



are zero at 




x 



0, TT . 



(13) 



Note that differentiating (12) gives: 




dx 



0 , x >(j+l)Ax or x < (j-l)Ax 

1/Ax , (j-l)Ax _< x jAx 
-1/Ax , jAx £ x £ (j+l)Ax 



(14) 



The left hand side of (13) is easily evaluated since only 3 terras in the 

sum are different from zero: 

TT 




dcp. dm. u. -Ax - 2u,Ax + u. ,-Ax 

_± _LL dx = i+1 

dx dx .2 

Ax 



(15) 



The right hand integral in (13) may be evaluated if f (x) is approximated 
in terms of the basis functions: 



f (x) 



N 



E 




y 



so that the integral becomes 




(i+1) Ax 



f j / cp^jdx . 

(i-*l)Ax 



( 16 ) 
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If £ = x - lAx is introduced the integral can be expanded into three in- 




cp^f (x)dx 



1 “ 1 J Ax 2 1 j Ax 2 1+1 I Ax 2 

-Ax -Ax 



(17) 

When these terms have been evaluated, (17) and (15) can be substituted 
into (13) which gives 



u (11 - 2u. + u. - 
l+l i i-I 



f _ + 4f . + f. - 
1+1 1 1-1 



Ax 



(18) 



This equation applies for 2 _< i <_ N-l and the equations for i = 1 and 
i - N are obtained by removing any terms in i = 0 or i = N+l . Equation 
(18) may be solved by Gaussian elimination 

Since each coefficient in this finite element expansion represents the 
solution at a certain point in space, it is convenient to compare (18) 
with finite difference forms of (6). The centered difference form (4) 
of this equation is 



U i<-1 ~ 2u t + u i-l 

Ax 2 



- £ 1 > 



(19) 



where u_^ = u(iAx) . The finite element equation (18) and the finite 
difference equation (19) are the same, except that the forcing term in 

(18) appears in a weighted average. When these equations are solved with 
a f(x), which is sinusoidal, the finite element form is considerably more 
accurate for the shorter wavelengths. 

In this example it appears that the spectral method is superior because 
the solution error is actually orthogonal to the basis functions. This is 
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not generally true with the finite element method because 



L(u^) depends on u^ ^ , u^ and • Each increase in 

N will normally change all of the solutions u_^ , whereas with the spectral 
method the original N amplitudes are not changed because they are already 
exact. However if the variation of f should require fine resolution in 
only a small area, the finite element method can easily be applied by letting 
Ax vary. In this case the spectral method would require more elements be- 
cause its spatial resolution is uniform. It is also clear that the finite 
element method can be used to design better finite difference equations. 

3 Time Dependence 

In the previous sections the Galerkin procedure has been applied to one- 
dimensional equations which are independent of time. The treatment of time 
variation is important for most meteorological prediction problems. Consider 
the following simplified equation: 




X (u) 



= 0 , 



( 20 ) 



where the operator Jtf may be nonlinear. Expand u(x,t) into a series as 
follows : 



u(x,t) 



N 

^ ' u (t)cp:(x) , 

j=i 3 J 



(= 21 ) 



where the coefficients u^(t) are functions of time and the basis functions 
C(X (x) are functions of x . Usually the Galerkin procedure is not applied 
to the time dependence because it is more convenient to ue finite differences 
in time. 

The Galerkin form of (20) is obtained by substituting (21) into 
(20), multiplying by qx(x) and integrating over the domain as follows: 
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N 



J c Pi C P i dx + / cp i( £( s u cp )dx = 0 , i=l,...,N . 

j = l a J J J 



( 22 ) 



This process gives N coupled ordinary differential equations in the co- 
efficients u j(t) • This set can be solved by introducing finite differ- 
ences in time. 

The importance of energy conserving finite difference schemes is well 
known. The Galerkin method leads naturally to energy conservation in equa- 
tions with quadratic energy variants. To show this, multiply (20) by u 
and integrate with respect to x : 



9(u 2 /2) 

3t 



dx = - 



u£(u)dx . 



(23) 



For an energy conserving system, the operator must satisfy the condition 

b 

u£(u)dx = 0 , (24) 



where u is any function which satisfies the boundary conditions. In this 



case (23) becomes 



/ u 2 /2 dx = 0 , 



(25) 



which shows the energy conservation for the exact equation. To demonstrate 
that the same result holds for the finite sum (21) , multiply the ith 
equation of (22) by u^ and sum from i=l to i=N : 



b N 

( 



N 



b N 

^ u i cp i ) a7 ( ^ “iV* 1 * = " f u i c f ) i ) / <Xj u i c P i 

i=l dt j-1 J 2 / 1=1 j=l 2 2 

a 



) . ( 26 ) 
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The integral on the right vanishes from ( 24; since the function given by 

21) satifies the boundary conditions. Therefore (26) can be written 



b 



N 




(27) 



a 



which expresses the energy conservation for the Galerkin approximation to the 
spatial variation. As with finite difference equations the actual degree of 
energy conservation will depend on the time differencing which is used in 



4 Barotropic Vorticity Equation with Fourier Basis Functions 

In this section the spectral method will be applied to the barotropic 
vorticity equation on the beta plane. Fourier basis functions are appropri- 
ate for the beta plane when the fields are periodic in x and y . The 
development of this section closely follows Lorenz (1960). The barotropic 
vorticity equation may be written: 



where \[j is the streamfunction. Suppose that the fields are periodic in 
both x and y so that 



With the beta plane geometry and the periodicity condition, the appropriate 
orthogonal basis functions are of the form: 



( 22 ). 



v 2 \p + k X ViJj • V(V 2 ^) + 6 dip/dx = 0 , 



(28) 



^(x + 2TT/kjy + 2ir/Z,t) = 4>(x,y,t) . 



(29) 



cp ( x >y) = e 

mn 



i(mkx+n£y) 



(30) 



These functions are eigensolutions of the equation: 



V cp + bcp 3 0 



(31) 
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where the eigenvalues are given by 



2 2 2 2 
b = (m k + n Jl ) 



(32) 



The streamfunction can be expanded in terms of this basis functions as 
follows : 



iKx,y,t) = 




m n 



c ( t ) e 
mn 



i (mkx+n£y) 



In order for if; to be real the coefficients must satisfy the condition 



C = C 
mn -m-n 



where ( ) indicates the complex conjugation. This can be shown by con- 

sidering only the m,n and -ra,-n . It is convenient to introduce the wave 
number vector M = mki + nilj and the radius vector R = xi + yj . The 
expansion for can now be written 



iM*R 



iKx,y,t) = ^ , C+(t) e Ari R . (33) 

M 

With the use of (31) and (33) the vorticity can be written 



v 2 <|> 



E -v -v 

(M-M) C+(t) 



(34) 



M 



The quantities which are required in the nonlinear term in (28) may be 
written: 

x \ -> 

•R 



V<J> 



V(V 2 ^) 



'ST' r iH*i 

Z—J lH C H 6 



H 






iL(L*L) C+ e 



iL*R 



(35) 
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The wavenumber vectors H and L are introduced because the sums must be 



multiplied together and rearranged. 

Now substitute the various sums [( 33), (3*0 and (35)] into (28) 

which gives: 



Y. C- e iL ‘ R +y y (L-L)k.HxL C+C+ e l(H+L) * 



R 



L H 



£ 



+ imBV , C+ e iL * R = 0 . 



C 36) 



The Galerkin method for this equation is similar to the method used in 
(22), except that the equation must be multiplied by the complex conjugate 
of the basis function since the basis function is complex. To carry out this 

*-+’ -f 

process multiply ( 36) by e and integrate over the periodic domain 

as follows: 

27r/k 2 tt / Si 



1 1 






(L-L)k-HxL CjjCj e 1 ( H+I *~M) • R | dydx = 0 > 



L H 



for each M in the original sum (33) . Each integral of the exponential 
function will vanish except when the exponent is zero. This leads to the 
following equation for each M t 



d<g taWiGJ 
dt 



" eo u , y 

M*M 

H. 



(M-H)» (M-H)k»HxM 

- 4 - 

M*M 



^-H C H ° * 



(37) 
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In the first two terms the contribution occurs for L = M and in the last 



terra for L = M - H . 

Equation (3) represents N ordinary differential equations, where 
N is the number of terms in the sum (33). The last term in the equation 
gives the interaction between different waves which comes from the nonlinear 
advection terra in (28). In particular wave M is affected by the inter- 

*f -f -f r « 

action of waves H and M - H . V/ hen the last term is dropped, (3) 
becomes a set of linear, uncoupled equations which can be solved to give 
the Rossby wave solutton. 

In section (3) it was pointed out that the Galerkin procedure pre- 
serves energy type univariants which arise from quadratic nonlinearities in 
the original equations. Equation (28) conserves both kinetic energy and 

mean square vorticity or enstrophy. The kinetic energy for the region can 
be written: 



27T/k 2-n/l 27t/k 2ir/fc 




where the product was obtained from (35) with different summations. 

*>■ 

The integral on the right is nonzero only when H = -M so that the energy 
can be written 




11 

M 




(38) 



where the condition C ^ = C ^ has been used in the last step. 

-M M 

The energy form in (38) is conserved (dK/dt = 0) by both the original 
vorticity equation (28) and the spectral form (37). The conservation for 
( 28 ) is easily demonstrated, and the conservation for ( 38 ) follows 
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from the development in section (3). An equation for the rate of change 
of energy in wave M can be obtained by differentiating with respect 

to t and by using (37). The resulting equation shows that the energy 
in wave M changes in proportion C+ times the amplitudes of pairs of 
interacting waves. Thus if Cg is maintained at zero, the energy flow out of 
the other waves to it must be zero. This shows in another way that energy 
will be conserved in any set of waves that might be selected for sum (33). 

S ince interactions outside of this set are neglected, aliasing cannot occur 
in a spectral model. This automatically eliminates the nonlinear computa- 
tional instability which occurs with finite difference equations. 

The set of ordinary differential equations (37) can be integrated 
numerically with one of the standard schemes. In fact Baer and 
Platzman (1961) noted that the linear terms in (37) can be treated exactly 
so that the only time differencing errors comes from the nonlinear terms. 

It is clear that the spectral method is much more accurate than most 
finite difference methods for the same number of degrees of freedom. In par- 
ticular, linear advecti.on that was examined is treated exactly 
by the spectral method provided that the initial field in resolved. Finite 
difference methods experience false dispersion since the short waves move 
too slowly. The spectral method has no aliasing because interactions in- 
volving shorter waves outside of the truncated set are excluded. On the other 
hand, the finite differencing falsely reflects interactions with shorter waves 
back onto longer waves. With the Arakawa Jacobian finite difference forms, 
this aliasing does not produce spurious energy, but it 

does cause phase errors in the interacting waves. In spectral models the 
most important error involves the neglect of Interactions with wave components 
which are outside of the original set. The neglect of these interactions 
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causes an error in the waves which are represented by the basis functions. 

Thus although the error in the original equation is orthogonal to the basis 
functions, the error in the solution will occur in the scales described by 
the basis functions. 

When the spectral method is applied to a vorticity equation such as 
(28), a Poisson equation for 9^/St does not have to be solved since the 
basis functions are eigensolutions of (31) . The Poisson equation must be 
solved at each time step with finite difference methods. The biggest draw- 
back to this form of the spectral equations is in calculating the nonlinear 
term which appears as the sura in (37). The coefficient preceding 

is called the interaction coefficient and it is usually computed 

M-H H 

just once and stored for use during the integration of the equation. The 

problem is that if there are N degrees of freedom the number of operations 

2 

needed to compute the nonlinear terms goes as N for this spectral model 
as compared with N for most finite difference methods. Thus for high 
resolution (large N), this form of the spectral method requires relatively 
larger computer time than finite difference methods. In a later section a 
method which avoids this problem will be presented. However the present 
method is very convenient for low-order models. Lorenz (1960) obtained some 
very interesting nonlinear solutions with a 3-component system. It can be 
seen from ( 37) that at least 3 waves are required for nonlinear interaction. 
5 Barotropic Vorticity Equation with Spherical Harmonics 

In this section the spectral equations will be formulated for barotropic 
motion on the sphere. The barotropic vorticity equation in spherical co- 
ordinates can be written: 

8 2. 1 rdip dV 2 ip dip 3V 2 iK 2Q ( 39) 

3t ^ d\ 3y J 2 3X U ’ V ' 

a a 
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where 



v 2 = itaV) 



(40) 



1-y 3A' 



In these equations X is the longitude and y = sincp , where cp is the 
latitude. The spectral method was first applied in spherical coordinates 
by Silberraan (1954) and the development of this section follows Platzman 
(1960). 

The appropriate orthogonal basis functions are 



Y (y, A) = P (y) e lmX , 
m,n m,n 



(41) 



where P 

m,n 

are defined 



denotes associated Legendre functions of the first kind which 
by 



p m ^ 

m,n (n-Hn; ! „n , 

z n : 



2.m/2 ,n+m _ 

d ,2 , N n 

(y -1) 



dy 



n-hn 



C42) 



These basis functions are spherical harmonics which satisfy the equation 



V 2 Y + b Y = 0 , 
m,n m,n 



where, the eigenvalues are given by 

b = n(n+l)/a^ . 



(43) 



(44) 



Here |m| is the planetary wave number and n-|m| is the number of zeros 
between the poles. Also n must be greater than or equal to |m| . These 
basis functions are orthogonal so that 




1 for (m',n') = (m,n) 
0 for (m’,n') 4 (m,n) 



(45) 
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The streamf unction can be expanded as follows: 



M |mj+J 

<Kx,y,t) =a 2 ) X'' ip (t) Y (X,y) 

/ j / i m,n m,n ’ 

m=-M n=|m| 



( 46 ) 



Since ip must be real ip must satisfy 

m,n J 



\ p = (-D ra ip* 

-m,n m, n 



(47) 



This condition was ’derived with the use of the relation P = (-l) ra P 

-ra,n m,n 

The coefficients ip can be obtained from the inverse transform: 

m, n 

27 T 1 



^n,m^ t ^ 47Ta1 



tp(A,y,t) Y dpdX 
m, n 



(48) 



The vorticity has the following expansion: 

M |m|+J 



V 2 ip 



E- 

ra=-M n=(m| 



(n+l) ip (t) Y 

m,n ra,n 



(49) 



which follows from (43) and (44). 

_The Galerkin method is applied by substituting (48) and (49) into 

* 

(39), multiplying by ^ and integrating with respect to p and A . 

When the conditions (45) are employed this equation reduces to: 



dip 



20m i 



dt 



n(n+l) ^ra,n n(n+l) r m,n 



(50) 



The nonlinear terms F may be written 

m,n 



M |m|+J 



M 



|m |+J 



m,n 



Z , Z X) , L <m,n;VWn 2 ) 

»]— M nf | m 1 1 m 2 =-M n 2 =|m 2 | 



( 51 ) 
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where the interaction coefficients are given by 






dP 



jtnjCnj+D-njCnj+D] | 



m 2* n 2 



dP 



-1 



- m«P 



TTl^ , ft 



2 m 2 ,n 2 dy 



]dy for m = m^+ tn 2 



= 0 



for m ^ m^+ m 2 



\ 




(52) 



In obtaining this result the subscripts 1 and 2 were used for expansions 

(48) and (49), respectively. This form for the interaction coefficients 

comes from the fact that F changes sign when ip and ^ are 

m,n m l ,n i m 2 ,n 2 

interchanged. 

Equation (50) has the same form as the prediction equation (37) for 
the Fourier basis function. However the spherical coordinate equation has 
more complicated interaction coefficients because of the integral involving 
the Legendre functions. It can be shown by the same method as before that 
energy is conserved, and Platzman (1960) has also shown that mean square 
vorticity is conserved. The spectral method applied to spherical (global) 
prediction has the advantage that there are no singular points whereas singular 
points often cause problems with finite difference models. The only major 
disadvantage in solving C^O) is in the large number of terms which come from 
the nonlinear terms. This problem will be treated in the next section. 
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6 . 



Transform Method 



In this section a new method for handling the nonlinear terms in (50) 
will be presented which avoids the use of interaction coefficients (see (51) 
and (52)) . This method was formulated independently by Orszag (1970) and 
Eliasen, Machenhauer and Rasmussen (1970), and it has been reviewed by Bourke, 
McAvaney, Puri and Thurling (1977). The problem with the interaction coeffi- 
cient method for computing nonlinear terms is that it requires multiplication 
of two series (together) which is very time consuming. The transform method 
sums the series at certain spatial grid points and these fields are multiplied 
together at each point to form the nonlinear terms. Then the nonlinear terms 
must be transformed back to spectral space. The usefulness of this process 
is enhanced by the existence of efficient transform methods. In spherical co- 
ordinates the fast Fourier transform is used in longitude and the Legendre inte- 
grals in latitude are evaluated by Gaussian quadrature. This method is far 
superior to the interaction coefficient method for the sphere. 

The nonlinear terms which must be transformed may be rewritten as follows: 



f<p,X) - ^llrr ^ - If - ^twClrr 72 *> - £<4t ? ! <»! • (53) 



2 l 3y 3A dX 3y 



2 L 3A 3y 



3y 3A 



It is now convenient to define the following quantities which are the X and 
cp velocity components multiplied by coscp: 



it - ~ CQS CD jW v rl 

U ~ a 3y ’ " a 3A 



(54) 



When these velocities are introduced into (53) it can be written as follows: 



F(y,A) = - -j[ — — 

a i - y 



r^ (w2 *> + i (vv2 *» • 



(55) 
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The velocity components (.54) can be computed from (46) at longitude- 



latitude grid points, and the vorticity can be obtained at the same points 
using 49), The details of the process will be given later. The products 
UV 2 ip and VV 2 i p can be calculated at each grid point and the resulting prod- 
ucts can be Fourier analyzed in X to give the following relations: 



m=M 

UV 2 ifJ = a V"' A (y) e imA , 

/ j m 
m=-M 

m=M 

VV 2 ^ = a B m (y) e imX . 

m=-M 



(56) 



The transform of F(y,X) is given by 



F 

m,n 




-imX 

e 



P F(y ,X) dydX . 
m, n 



(57) 



The X integration in (57) can be carried out by substituting (36) into 
(55) and by inserting the result in (57), which gives: 



F 

m,n 



1 

2 





i-y 



2 



A P 
m m,n 




p 

m,n 



] dy • 



The second term can be integrated by parts which gives 



F 

m,n 



1 

2 




-1 



r im 



i-y 



A P 
m m,n 



dP 



m,n 



m dy 



] dy , 



(58) 
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where the condition — 0 at y = ± 1 was used to simplify the integral. 

This condition follows since V is equal to the actual velocity times coscp. 

The form of F ^ given by C.58) is superior to the earlier form because 

only the known function P is differentiated. 

ra,n 

The integrand in (58) is a polynomial in y and the integral can be 
evaluated following Eliasen et al. (1970) by the Gaussian quadrature formula. 
If the integrand is denoted by Q(y) , the formula gives the following 
expression for F 



m,n 



R 



f = \ g. 

m,n 2 / a 



r K) «<V> ■ 



159 ) 



k=l 



In (59) the sununation is carried over K values of y , where the y 's 

K. K 

(K) 

are roots of the Legendre polynomial P and ' are the corresponding 

W, J\ K. 

Gauss coefficients. The formula is exact for any polynomial of degree smaller 
than or equal to 2K-1 (see ). Thus apart from 

roundoff errors, no approximation is introduced by computing the integral when 
a sufficiently high value of K is used. The maximum degree of Q(y) can be 
most easily obtained from (52) . 

Before discussing this process for treating the nonlinear terms in more 

detail,~it is necessary to determine the relation between J and M which must 

be defined in the sum (46). In rhomboidal truncation J = M , so that each 

latitudinal mode has the same numbej: of waves in longitude. With triangular 

truncation J = M - |m| so all basis functions which have the same scale 
2 

n(n+l)/a , are either retained or dropped. Thus the mode with the smallest 
latitudinal scale has the largest longitudinal scale. Most meteorological 
models use the rhomboidal truncation in part because it gives better longitudinal 
resolution. In the remainder of this development, the rhomboidal truncation 
will be used. 
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In order to construct the fields (56) it is necessary to obtain U 



and V from ip . First expand U and V into these sums: 



M | ml+M+l 

u = ±y V u y , 

a m,n in , n 

m=-M n= jm| 



M |m|+M 



”iE E. 

m=-M n= I rn I 



V Y 
m,n m,n 



(60) 



The following relations will be useful in evaluating (54) : 

? 9Y 

(U -1) = n D - Y - (n+1) D Y 

9y m,n+l m,n+l m,n m,n-l 



(61) 



8Y 



ax 



= in Y 



m,n 



2 2 2 H 

where D = [(n -m )/(4n -1)]^ . The final expressions for U and 
m,n m,n 

^ can be obtained by substituting (46) and (60) into (54), using 
(61) and by applying the orthogonality condition (45) : 



U - (n-1) D ip - (n+2) D ip _ , 

m,n m,n m,n-l m,n+l m,n+l 



V = im il; 
m,n ra,n 



(62) 



Note that the expansion for U as given in (60) must extend one degree 



above that defined for ip , since nonzero values of U 
implied by nonzero values of ^| m | | m |+M * 



m| , |m|+M+l 



are 



The quantities U , V and V 2 ip can now be evaluated at points 



= 2ttJ/N , = arcsin 



28 



where j = 1 , . . . ,N and k = 1,. 
latitudes. Consider for example 



V( YV 




. .,K . The cp^'s are called the Gaussian 

V(A.,y ) which can be written: 

J k 

|ml+M ~1 



E lm \\) P (y ) 
m,n m,n k 



(63) 



n= | m 



with the use of (41), (46) and (62). Similar expressions can be 

written for U(A.,y ) and V 2 i^(A.,y ) . The outer summation can be carried 
J k J k 

out very efficiently with the use of the fast Fourier transform method which 

was developed by Cooley and Tukey (1965) . The number of operations required 

for the fast Fourier method applied over N points is of order N N , 

2 

while for the direct method order N operations are required. The fast 

Fourier transform method is clearly much faster than the direct method for 

larger values of N. The next step is to compute UV 2 ip and W 2 ip at each 

grid point. After these products have been computed, the Fourier transforms 

must be calculated to give A and B for use in (56). For example, 

mm 

using the discrete Fourier transform: 

w=sE e 3 <uv 2 w j k - (64> 

j-i 

where -M _< m _< M . A similar expression is obtained for B^y^) . The 
fast Fourier transform can also be used here to save time. 

It is important to choose N large enough to avoid aliasing when the 
products are transformed back to wave number space as in equation (64) . 
Orszag (1969), (1970) suggested that N = 4M would be needed, but later 
Orszag (1971) and Machehauer and Rasmussen (1972) showed that N = 3M+1 was 
adequate to provide alias-free transforms. 
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can be computed exactly 



Now that A (u. ) and B (u, ) are known, F 

ra k ra k m,n 

from (59) if the degree of the polynomials is less than or equal to 2K-1. 

The maximum degree can be determined from (52) by noting that n is a 

polynomial of degree n and by considering these selection rules for the 

interactions: m = m^ + m^ , \Si^ - H^\<Z < • T ^ e conclusion which 

is given in Bourke et al. (1977) is that the maximum degree is 5M-1, so that 

the number of Gaussian latitudes is K > 5 M/2 . 

This method of computing F is more efficient than the interaction 

m,n 

coefficient method and it requires much less computer storage. The number of 

calculations required for the interaction coefficient method is of order (M^) 

3 

while for the transform method it is of order (25 M ) [see Bourke et al. 
(1977)]. It will be shown in the next section that the transform method is 
more efficient for even a moderate value of M and this advantage increases 
rapidly with M . 

7 Spectral Model of Shallow Water Equations 

In this section the spectral method will be extended to the primitive 
equations and it will be demonstrated that semi- implicit differencing can be 
applied with little extra effort. The shallow water equations in spherical 
coordinates will be used to demonstrate the procedure following Eliasen et al. 
(1970)" and Bourke (1972). The equation of motion and the continuity equation 
can be written: 

3v v«v 

= - (C+f) k x V - V(4>* + , (65) 

f£'- " V*<f>'V - <f>6 . (66) 

This form of the equation of motion will simplify the derivation of the vor- 
ticity and divergence equations. Note that the geopotential has been split 

_ i 

into a mean $ » and a departure <f> , which will facilitate the implementa- 

tion of semi- implicit time differencing. 
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The velocity is broken into rotational and divergent parts as follows: 



V = k x Vi|> + Vy = (U/coscp) i + (V/coscp) j . (67) 

The modified components U and V will also be used here. Mow form the vor— 
ticity and divergence equations by taking V* and k*Vx of (65) which 
gives : 



|f = - v-(c+f) V , 



( 68 ) 



|f = k-Vxfo+f) V]- V 2 (<f> * + ^) . 



(69) 



The vorticity and divergence become 
C » v 2 ^ , 6 = V 2 x . 



(70) 



In spectral models it is convenient to replace the equation of motion by 
the vorticity and divergence equations because the relations (70) are sim- 
plified when spherical harmonics are used as basis functions. This form of 
the equations is also more convenient for application of semi-implicit dif- 
ferencing. 

The vorticity equation (68) and the divergence equation (69) can 
now be expanded with the use of (67) and (70) to give: 

A v 2 ^ [A (u v 2 i|>) + coscp^ (vv 2 ip) ] 

a cos cp ^ 

- 2fl (sincp V 2 X + V/a) , (71) 

£ y2 X - —V (W 2 *) - cos„|- (UV 2 *)] 
a cos cp ^ 

2 2 

+ 2J2( sincp V 2 ip - U/a) - V 2 ( U + | + 4> * ) . (72) 

2 cos cp 
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Similarly the continuity equation (66) becomes: 



9 £' 

9t 



1 

2 

a cos cp 




W') + cos cp ^ 



The two components of (67) can be written: 



_ _ coscp 9 iJj 9x. 

a 9cp a 9A * 



v = 1 M + cos cp 

a 9X a 9cp 



( 73 ) 



(74) 

(75) 



Equations (71), (72) and (73) are the predictive equations for , x 

and <p ' and (.74) and (75) are diagnostic expressions for U and V . 
The nonlinear terms in these equations are in a convenient form for the 
transform method which was presented in the last section, since the multi- 
plication can be performed at the grid points before differentiation. 

Each of the dependent variables is expanded in terms of the spherical 
harmonic basis functions as follows: 




( 78) 

These expansions are for the rhomboidal wave number truncation. Equations 
(74) and (75) are transformed in the same manner as equations (54) were 
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in the last section and the result is 



U 

m,n 



(n-1) 



D ip . 
m,n m,n-l 



- (n+2) D ip + irav 

m,n4-l r m,n+l A m,n 



V 

m,n 



=-(n-l)D x i +(n+2)D .y . + imili 

nijn^m^n-l m,n+l A m,n+l r ra,n 



(79) 



Note that the expansions for U and V must extend one degree above the 
expansions for ip and x • 

The quantities needed for the nonlinear terms are obtained by evaluat- 
ing the sums in (76), (77) and (78) at equally spaced points in longi- 

tude and at Gaussian latitudes. The required products are computed at each 
point and the products are then Fourier transformed in longitude as follows: 



M 



M 



UV 2 if> 



= a 



A e 
ra 



imA 



W 2 ip 



= a 



B e 
m 



imA 



m=-M 



m=-M 



v *' - a3 £ c „ elmX - v +* d . 

m=-M m=-M 



imA 



. 2 . 2 



M 



2 / j m 



imX 



m=-M 



(80) 



(81) 



(82) 



The spectral equations are formed by substituting (76), (77), (78), 

(80), (81) and (82) i nto the system (71)- (73) and multiplying each 

'k 

equation by Y and integrating over the domain. With the use of the 

m,n 

orthogonality condition (45) the equations finally reduce to the following 
set: 



-n(n+l) 



8i|^m,n _ 
3t 




1-Ji 



TT [imA P 
2 m m,n 



- B 



dP 

B 

m dy 



m,n 



] d y 



+ 2fi[n(n-l)D Y n , + 
m , n m y n X 



(n+l) ( n+2 ) D m>n+1 X m>n+1 



- V ] 
m,n 



(83) 
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1 



-n(n+l) ^ - i 



/ dP 

-i-r-tiraB P + A ~^] dy-2ft[n(n-l)D J> 

m m,n d dp m,n m,n-l 

-j. 



+ (n+1) (n+2)D , .. + U ] + n(n+l)(E + <J> ) » 

m,n+l m,n+l m,n m,n m,n 



(84) 



if 1 

^■-2 TV 11 

^1 



dP 



[imC P - D . n m?n ] dy + (f>n(n+l) Y , (85) 

m m,n m dy m,n 



where 



m 



J m,n 2 



-1 



i-y 



0 p dy . 

1 m,n 



( 86 ) 



The integrals are evaluated by the Gaussian quadrature formula as before, 
but this time (5M+l)/2 Gaussian latitudes are required. As before the 
required number of longitudinal grid points is 3M+1. 

Bourke (1972) compared the efficiency of the transform method to the 
interaction coefficient method for this model. Figure 2 shows the computer 
time required per time step for the two methods as a function of the trunca- 
tion number M. The figure shows clearly that even for M = 15 the transform 
method is an order of magnitude faster than the interaction coefficient method. 

In fact the interaction coefficient method becomes almost intractable for M 

« 

much larger than 15. At M = 15 there are over 500,000 interaction coeffi- 
cients. 

The system (83)- (85) is very convenient for the application of semi- 

implicit time differencing. All terms are evaluated explicitly except that 

<J> in (84) and Y in (85) are treated implicitly. These two 
m,n ra,n 

equations are easily solved for <j> m ^(t+At), and equations (83) and ( 84) 
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Figure 



2 

-Computation time per time step (s) as a function of 
spectral resolution. Integrations of a global spectral model 
employing a transform method and employing the interaction 
coefficient method are compared. 



can then be solved explicitly. In contrast finite difference models require 
the solution of a Helmholtz equation for <f)(t+At) , at every time step 
Thus in spectral primitive equation models a much longer time step can be 
used with almost the same computational effort per time step. 

The introduction of the transform method and semi-implicit differencing 
have made the spectral primitive equation models competitive with finite 
difference models for global prediction. _ The procedures used in this section 
are easily extended to baroclinic models as has been done by Bourke et al. 
(1977), Machenhauer and Daley (1972) and Hoskins and Simmons (1975). Com- 
parisons have shown that as good or better forecasts can be made with global 
spectral models than with finite difference models which use the same amount 
of computer time (Doron et al. (1974) and Daley, Girard and Simmonds (1976). 

It should be pointed out that energy is not exactly conserved in this 
model even with continuous time variation. This is because the kinetic energy 
for the shallow water equations is proportional to <j>V*V which is a cubic 
energy form, and consequently the analysis of Section 3 does not apply. 
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However the nonlinear terms are computed very accurately in spectral models 
and experience shows that the energy is in fact very nearly conserved. 
Bourke (1972) integrated the model which was developed in the section for 
116 days, and obtained an energy change of only 2 percent. 

8 Advection Equation with Finite Elements 

In this section the finite element method with linear elements will be 
applied to the advection equation 



8u Su 

St C Sx 



= 0 . 



(87) 



This equation has heen treated extensively wfth various finite differs 

V S 

ence schemes. The Galerkin equation is obtained by setting JL - c — in 
( 22) which gives 




The linear basis functions cpj (x) are defined by (.12) and a typical one 
is shown in Fig. 1. In this application u is periodic so that the basis 
functions must satisfy cp Q = Cf^ and cp^ = qp^ . 

The first term in (88) can be evaluated from (17) which is of the 
same form, and the second term can be computed with the use of ( 14) . The 
resulting equation with i = m can be written: 



- du , - du du 

1 ( _nrfl 4 _m _ 

6 Mt dt dt 



1 * U 1 
m-1) + c j afl m-1 = 

2Ax 



(89) 



The advection term is the same as is obtained from centered differencing, 
but the time derivative appears as a weighted average over three points. It 
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will be seen later that this greatly increases the accuracy of the solution. 
Now apply leapfrog time differencing which gives the following equation: 



12At (u m+l,n+l U nH-l,n-l +4 (u m,a+l _U m,n-l )+U m-l,n+l U m-l,n-l )+ 2Ax^ U m+l > n'‘ U in-l ,n ) 



(90) 

The stability and phase error can be investigated by substituting u = 

m , n 

A exp [ i (yAxra + aAtn) ] into (90). This yields 

sinaAt = -(cAt/Ax)(3 sinyAx)/(2 + cosyAx) . (91) 

The solution is stable (i.e. (sinaAt) £ 1) if 

3cAt/Ax[sinyAx/(2 + cosyAx)] < 1 . 

max — 



The bracketed term is a maximum when yAx = 120° , so that the stability con- 
dition becomes 

|cAt/Ax| £ 1//3 . (92) 

This criterion is considerably more restrictive than the CFL condition which 

arises from the leapfrog finite difference scheme. However it will be shown 

that ' (90) gives even better phase speed than the fourth-order leapfrog 

scheme for which the computational stability criterion is jcAt/Ax| £ 0.73 

Thus it is not unreasonable that ‘the leapfrog finite element scheme would 

have a more restrictive computational stability criteria. 

The finite element formula with leapfrog time differencing is actually 

implicit, since the new value u cannot be obtained explicitly from the 

ra, n+1 

earlier time values. Thus it seems reasonable to use a fully implicit form 
which does not have the timestep restriction (92). Consider the following 
time difference approximation to (89): 
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C 93) 



(u . -u + 4(u ,-u ) + u _ . — — u , ) 

6At nrKL,n+l nrKL,n m,n+l ra,n m-l,n+l ra-l,n 



+ 7r( u ^.i ,,-u i . i + u , i “ u i ) = 0 . 

4 Ax nrt-l,n+l m-l,n+l m+l,n m-l,n 



This fully implicit scheme can be shown to be neutral for all timesteps, 
and it should require about the same effort per time step as (90). For 
this reason implicit time differencing schemes are often desirable when 
finite elements are used. 

The phase speed for the leapfrog scheme is given by 



Cp = - a/y 



1 

yAt 



arcsin 



c At 3 sinyAx . 
Ax 2+cosyAx 



(94) 



If At/Ax and y are fixed, this expression approaches c as At 0 , 
which shows that the solution converges. If At is small in comparison to 
Ax/c 9 this formula reduces to 

_ c 3 sinyAx _ __c sinyAx , 

C F = vtx 2+cosuta ‘ uAx u . 2/ 3 sln 2 (M4x/2)] ' 

Table 1 contains c/c from ( (95) for typical values of L . 

F 



L 


2Ax 


3Ax 


4Ax 


6 Ax 


FEM 


0 


0.83 


0.96 


0.99 


4th order 


0 


0.61 


0.85 


0.96 



Table 1: c/c^ for the FEM solution and for 4th order space 

differenced scheme for various wavelengths L. 

The table also includes the ratio for the fourth order scheme from the limit 
for small At . The finite element formula (95) can be expanded 
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in pAx which leads to an error that is of order (pAx)^* . Table 1 shows 



that although the linear finite element equation and the fourth order finite 

difference equation have the same order of truncation error, the finite 

element equation is much more accurate. At L = 3Ax the finite element 

solution gives only 17% error in phase speed, while the fourth-order finite 

difference gives 39%. However for L = 2Ax , c = 0 , which indicates that 

r 

the finite element computational group velocity is very large for this wave- 
length. This can be shown by differentiating as follows: 



When L = 2Ax (pAx = tt) this formula gives G = -3 c which is much larger 
than the -(5/3) c that occurs with fourth-order differencing. This suggests 
that small scale noise will propagate very rapidly in finite element models. 

This tendency toward noisiness has been observed in various finite element 
models. The degree of accuracy indicated above for the finite element model has 
been verified by Cullen (1973) in a two-dimensional advective problem. It 
should be noted that although the FEM gives a solution for all values of x 
in the range considered, the high accuracy is only obtained at the nodal 
points since the fields are assumed to be linear between nodal points. In 
the next section the method will be applied to the barotropic vorticity equa- 
tion. 




[2cospAx -f 1] 

2 * 

(cospAx + 2) 



C dp 



( 96) 
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9 . Barotropic Vorticity Equation with Finite Elements 

In this section the finite element method will be applied to the non* 
linear barotropic vorticity equation in a two-dimensional domain. The basis 
functions will be linear functions on triangular elements. The barotropic 
vorticity equation can be written 



where 



= -k x v^-vn , 



n = f(y) + v 2 ip , 



( 97) 

0 98) 



is the absolute vorticity. 

Following Fix (1975) both ip and p are expanded in terms of the basis 
functions cp^ (x,y) as given below: 



N 

tKx,y,t) = i^(t) cpj (x,y) , 

3=1 

N 

n(x,y,t) = ^^ n^(t) cpj(x.y) . 

3=1 



( 99) 



( 100 ) 



When the Galerkin method is applied to ( 98) the following is obtained: 

dA = ff 0j(t) j'j'cp i (x,y) cpj (x,y)dA , 

for i=l,..,N. Since linear basis functions will be used it is necessary to 
integrate the left hand side by parts, which gives: 

N r F 



Xi tjj J Vcp-VcpjdA =J J cp.f(y)dA - 22 n j I I ’ 



(• 101 ) 



for i=l,2,...,N. 



bO 



The boundary terms which arise from the integration by parts were set to zero, 
by assuming that either ip is periodic in space or that there is no flow 
normal to the boundaries (i. e. ,kxVi^*n = 0, where n is a unit vector normal to 
the boundary) . Now apply the Galerkin method to the vorticity equation ( 97) , 

which leads to the following form: 



^2 i = _ ^2 X) ,*//■ ^ * x Vcp j * Vc ^ dA » ( 102) 



for i=l,...,N. This equation is of the same form that was obtained with the spectral 
model, but the nonlinear term requires much less effort because the only cp ? s 
which interact are those which are physically adjacent. 

The equations ( 101) and ( 102) conserve both mean square vorticity 

(enstrophy) and kinetic energy. The enstrophy conservation can be shown by 
multiplying ( 102) by r)^ and summing over i . When the sumations are 

taken under the integrals, the form ( 26) is found. Since the integral of 

-f 2 

Dk x vanishes, the conservation of 0 / 2 follows directly. The 

kinetic energy change can be examined by first differentiating ( 101) and 

substituting the result into ( 102) which gives: 




( 103) 



for i=l,2,...,N. Multiply this equation by and sum over i . The 

resulting equation is again of the same form as ( 26) and the left hand side 
is the derivative of the total kinetic energy. Since the integral of 
ijdc x Vi|/*Vn is zero, the energy is conserved. These results are not dependent 
on the particular basis functions which are employed. 

The systems of equations ( 101) and ( 102) can be written in matrix 

forms which are more convenient for solution. Let and q be column 



vectors of the values of and > respectively. Then ( 101) takes 



the form: 

= Q , ( 104) 

where the elements of the matrix K are 

K = If Vcp. • V C p i dA , (105) 

and Q is a column vector of the right hand side of ( 101). Similarly, 

system ( 102) becomes 

M - J , ( 106) 

where the elements of M are 



M ±j = // cp i cp j dA , (107) 

and J is a column vector of the right hand side of ( 102). 

The solution procedure will be illustrated for the case where leapfrog 
time differencing is used in ( 105) which leads to the equation: 



MAn = 2At J , 
n 



( 108) 



where An - Vl - Vl 



The matrices K and M are computed initially 
and stored for later use. The equations can be integrated beginning with 



ip. , D. , and n. • The right hand side of ( 108) can be computed 

J,n j,n-l J,n 

from ip. and n. , and that equation can then be solved for An. . This 
j» n j» n M j 

increment can be added to r], - to obtain n. . With these values the 

J »n-l J,n+1 

right hand side of ( 104) can be computed, and ( 104) can be solved for 

, and the process can be continued. In this procedure it is necessary 
to invert the matrices K and M during each time step. These matrices are 
very sparse since only adjacent elements interact. In some cases direct 



methods can be used, but iterative methods are much more flexible. 

Cullen (1973) has shown that the two dimensional advective stability 
criterion for linear elements is 

C 109) 

where d is the distance between nodal points. This is consistent with the 
one-dimensional result (92), because the step from one to two dimensions 
is usually achieved by replacing the grid size with d//2 . In this applica- 
tion f c | would correspond to the maximum velocity in the domain. Since the 
condition (108) is rather restrictive for At and since two matrices must 
be inverted per time step it may be worth while to use a fully implicit form 
similar to (93) . 

The natural generalization of the tent function in one dimension to two 
dimensions is a basis function which is composed of triangular elements. On 
each triangle the function varies linearly from 0 at two vertices to 1 at the 
third which is the nodal point for the basis function. Figure 3 shows how 
a typical basis function qx is constructed on a rectangular grid of nodal 
points. This function is the sum of the six plane surfaces that are associated 
with each triangle. The basis functions can be equally well constructed when 
the nodal points are irregularly located, and it is not necessary to have six 
triangular elements in the construction. 

The elements in the matrix equations (104) and (106) are obtained by 
evaluating the integrals in equations (101) and (106). These integrals 
can be reduced to a series of integrals over triangles such as are shown in 
Figure 3. Within each triangle any point is affected by only the three 
basis functions which have nodal points at the three vertices of the triangle. 
Zienkiewicz (1971) and Desai and Abel (1972) describe a convenient procedure 
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for evaluating the integrals over each triangle. This involves introducing 
triangular coordinates which vary linearly across each triangle in the same 
manner as the basis functions. The integrals can then be evaluated quite 
generally. 

A. rigorous mathematical analysis of the finite element method is given 
in the book by Strang and Fix (1973) . The stability and convergence of the 
method are discussed in considerable detail. Most finite element applications 
are based on a variational formulation rather than the Galerkin approach which 
has been used here, although the Galerkin method is most appropriate when time 
dependence is included. Pinder and Gray (1977) developed the finite element 
metood with the Galerkin approach, and gave applications in hydrology which 
has similar equations to those which occur in numerical weather prediction. 

The finite element method has been applied to atmospheric prediction with 
the primitive equations in shallow water form. Cullen (1974) and Hinsman (1975) 



Fig. 



3. 




Construction of the basis function 
rectangular array of nodal points. 






on a 



carried out global forecasts with these equations using linear basis functions 
on triangles as discussed in this section. The elements were efficiently 
arranged so that the area of each element was almost the same over different 
parts of the globe. Most global finite difference models have a large varia- 
tion in grid size between the equator and the pole, and consequently are not 
very efficient. 
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Staniforth and Mitchell (1977) reformulated the shallow water equations 
in terms of the vorticity and divergence as was done in section 6.6 for the 
spectral model. In this form semi-implicit time differencing can be applied 
easily, which allows a much larger time step. This is very important since 
the finite element method generally requires more computer time per time step. 
Staniforth and Mitchell also found very little noise generation in their fore- 
casts, whereas many finite element primitive equation models tend to generate 
small scale noise if no smoothing is used [Cullen (197^)]. 

The finite element method when applied to meteorological equations gives 
very accurate phase propagation and also handles nonlinearities very well. 

The main drawback to the use of the method is the requirement that an equation 
solver must be applied to invert a large matrix at every time step for every 
variable. The development of flexible exact solvers for these matrices is of 
great importance. The finite element method can easily be applied to variable 
resolution problems, but some finite element models do tend to produce noise 
probably as a result of the large spurious group velocity for the shortest 
wave. However, the formulation of Staniforth and Mitchell (1977) seems to 
reduce this problem considerably. Schoenstadt (1978) has shown that noise is 
generated in finite element models where all variables are carried at the same 
models points. When the variables are staggered at different model points:or 
when the vorticity and divergence equation are used this problem can be 
avoided. The general procedure used by Staniforth and Mitchell (1977) appears 
to be superior because semi- implicit differencing can be easily iraplimented, 
and the forecasts are not noisy. 
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